# Josh Fjelstul
# Cliff Carrubba
# Emory University

# libraries
library(ggplot2)
library(gridExtra)

# clean workspace
rm(list = ls())

# clear console
cat("\014")

#################################################
# read in data
#################################################

# set working directory

# read in data
dat <- read.csv("Fjelstul_Carrubba_APSR_replication_data.csv", stringsAsFactors = FALSE)

#################################################
# plot theme
#################################################

plot_theme <- function (font = "Helvetica", size = 14) {
  base_family <- font
  base_size <- size
  half_line <- base_size/2
  theme(line = element_line(colour = "black", size = 0.5, linetype = 1, lineend = "butt"),
        rect = element_rect(fill = "white", colour = "black", size = 0.5, linetype = 1),
        text = element_text(family = base_family, face = "plain", colour = "black", size = base_size, lineheight = 0.9, hjust = 0.5, vjust = 0.5, angle = 0, margin = margin(), debug = FALSE),

        axis.line = element_line(colour = "black"),
        axis.text = element_text(size = rel(0.8), colour = "black"),
        axis.text.x = element_text(margin = margin(t = 8), vjust = 1),
        axis.text.y = element_text(margin = margin(r = 8), hjust = 1),
        axis.ticks = element_line(colour = "black"),
        axis.ticks.length = unit(8, "pt"),
        axis.title.x = element_text(margin = margin(t = 16, b = 8)),
        axis.title.y = element_text(angle = 90, margin = margin(r = 16, l = 8)),

        legend.background = element_rect(colour = NA),
        legend.key = element_rect(fill = "grey95", colour = "white"),
        legend.key.size = unit(1.2, "lines"),
        legend.key.height = NULL,
        legend.key.width = NULL,
        legend.text = element_text(size = rel(0.8)),
        legend.text.align = NULL,
        legend.title = element_text(hjust = 0),
        legend.title.align = NULL,
        legend.position = "right",
        legend.direction = NULL,
        legend.justification = "center",
        legend.box = NULL,

        panel.background = element_blank(),
        panel.border = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.ontop = FALSE,

        strip.background = element_blank(),
        strip.text = element_text(colour = "grey10", size = rel(0.8)),
        strip.text.x = element_text(margin = margin(t = half_line, b = half_line)),
        strip.text.y = element_text(angle = -90, margin = margin(l = half_line, r = half_line)),
        strip.switch.pad.grid = unit(0.1, "cm"),
        strip.switch.pad.wrap = unit(0.1, "cm"),

        plot.background = element_rect(colour = "white"),
        plot.title = element_text(size = rel(1.2), margin = margin(b = half_line * 1.2)), plot.margin = margin(half_line, half_line, half_line, half_line),
        complete = TRUE)
}

#################################################
# Figure 2: Results for Hypothesis 3
#################################################

# model with effectiveness, support, and influence
b <- c(.0082373, .0117786, .0155704, .0196099, .0238924, .02841, .0331523, .0381059, .0432545, .0485786, .054056, .0596615, .0653671, .0711425, .0769552, .0827707)
se <- c(.0070539, .0065043, .0060044, .0056146, .005413, .0054782, .0058576, .0065482, .0075089, .0086873, .0100365, .0115191, .0131058, .0147732, .0165016, .0182737)
var <- seq(-1.25, 2.5, 0.25)

# plot data
plot.dat <- data.frame(b = b, se = se, var = var)
plot.dat$ub <- plot.dat$b + 1.96 * plot.dat$se
plot.dat$lb <- plot.dat$b - 1.96 * plot.dat$se

# make plot
plot1 <- ggplot() +
  geom_ribbon(data = plot.dat, aes(x = var, ymin = lb, ymax = ub), alpha = 1, fill = "gray75") +
  geom_hline(yintercept = 0, color = "black", lwd = 0.5) +
  geom_line(data = plot.dat, aes(y = b, x = var), color = "black") +
  plot_theme(size = 14) +
  ylab("Marginal Effect of Support") +
  xlab(NULL) +
  scale_x_continuous(limits = c(-1.25, 2.5), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-0.02, 0.12), expand = c(0, 0)) +
  theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        plot.title = element_text(hjust = 0, margin = margin(b = 20)))

# make plot
plot2 <- ggplot(data = dat, aes(influence)) +
  geom_density(color = "black", fill = "gray75", alpha = 1) +
  plot_theme(size = 14) +
  xlab("Influence") +
  ylab("Density") +
  scale_x_continuous(limits = c(-1.25, 2.5), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 1), breaks = c(0, 1), expand = c(0, 0)) +
  theme(plot.margin = unit(c(0.25, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 35, b = 0, l = 0)))

# model with effectiveness, distance, and influence
b <- c(0.0025394, -0.0011272, -0.0051111, -0.0094137, -0.0140332, -0.0189648, -0.0242, -0.0297264, -0.0355279, -0.0415842, -0.047871, -0.0543597, -0.0610181, -0.0678103, -0.074697, -0.0816358)
se <- c(0.0051879, 0.0046617, 0.0042795, 0.0041668, 0.0044334, 0.0051028, 0.0061158, 0.0073927, 0.0088698, 0.0105019, 0.0122566, 0.0141087, 0.0160362, 0.0180192, 0.0200382, 0.0220742)
var <- seq(-1.25, 2.5, 0.25)

# plot data
plot.dat <- data.frame(b = b, se = se, var = var)
plot.dat$ub <- plot.dat$b + 1.96 * plot.dat$se
plot.dat$lb <- plot.dat$b - 1.96 * plot.dat$se

# make plot
plot3 <- ggplot() +
  geom_ribbon(data = plot.dat, aes(x = var, ymin = lb, ymax = ub), alpha = 1, fill = "gray75") +
  geom_line(data = plot.dat, aes(y = b, x = var), color = "black") +
  geom_hline(yintercept = 0, color = "black", lwd = 0.5) +
  plot_theme(size = 14) +
  ylab("Marginal Effect of Distance") +
  xlab(NULL) +
  scale_x_continuous(limits = c(-1.25, 2.5), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-0.14, 0.02), expand = c(0, 0)) +
  theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        plot.title = element_text(hjust = 0, margin = margin(b = 20)))

# make plot
plot4 <- ggplot(data = dat, aes(influence)) +
  geom_density(color = "black", fill = "gray75", alpha = 1) +
  plot_theme(size = 14) +
  xlab("Influence") +
  ylab("Density") +
  scale_x_continuous(limits = c(-1.25, 2.5), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 1), breaks = c(0, 1), expand = c(0, 0)) +
  theme(plot.margin = unit(c(0.25, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 42, b = 0, l = 0)))

# save plot
pdf("Fjelstul_Carrubba_APSR_figure2.pdf", width = 9, height = 7)
grid.arrange(plot1, plot3, plot2, plot4, ncol = 2, nrow = 2, widths = c(4, 4), heights = c(4, 2))
dev.off()

#################################################
# Figure 3: Results for Inference 1
#################################################

# model with effectiveness, support, and influence
b <- c(-0.0411386, -0.0363879, -0.0313934, -0.0261524, -0.0206628, -0.0149231, -0.0089321, -0.0026893, 0.0038052, 0.0105506, 0.0175457, 0.0247884, 0.0322759, 0.0400047, 0.0479708, 0.056169, 0.0645935, 0.073238, 0.0820948, 0.091156, 0.1004125, 0.1098544)
se <- c(0.0118195, 0.0111524, 0.0104758, 0.0097903, 0.009098, 0.0084027, 0.0077115, 0.0070357, 0.0063938, 0.005814, 0.0053376, 0.005019, 0.0049164, 0.0050716, 0.0054901, 0.0061435, 0.0069884, 0.0079836, 0.0090963, 0.0103022, 0.0115832, 0.0129258)
var <- seq(-3, 2.25, 0.25)

# plot data
plot.dat <- data.frame(b = b, se = se, var = var)
plot.dat$ub <- plot.dat$b + 1.96 * plot.dat$se
plot.dat$lb <- plot.dat$b - 1.96 * plot.dat$se

# make plot
plot1 <- ggplot() +
  geom_ribbon(data = plot.dat, aes(x = var, ymin = lb, ymax = ub), alpha = 1, fill = "gray75") +
  geom_hline(yintercept = 0, color = "black", lwd = 0.5) +
  geom_line(data = plot.dat, aes(y = b, x = var), color = "black") +
  plot_theme(size = 14) +
  ylab("Marginal Effect of Effectiveness") +
  xlab(NULL) +
  scale_x_continuous(limits = c(-3, 2.25), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-0.08, 0.16), expand = c(0, 0)) +
  theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        plot.title = element_text(hjust = 0, margin = margin(b = 20)))

# make plot
plot2 <- ggplot(data = dat, aes(support)) +
  geom_density(color = "black", fill = "gray75", alpha = 1) +
  plot_theme(size = 14) +
  xlab("Support") +
  ylab("Density") +
  scale_x_continuous(limits = c(-3, 2.25), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 0.6), breaks = c(0, 0.6), expand = c(0, 0)) +
  theme(plot.margin = unit(c(0.25, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 33, b = 0, l = 0)))

# model with effectiveness, distance, and influence
b <- c(0.0854415, 0.0758935, 0.06652, 0.0573237, 0.0483067, 0.0394711, 0.0308186, 0.0223508, 0.0140691, 0.0059746, -0.0019318, -0.0096495, -0.0171778, -0.0245166, -0.0316656, -0.0386251, -0.0453952, -0.0519764)
se <- c(0.0110932, 0.0096579, 0.0083142, 0.007086, 0.006011, 0.0051482, 0.0045776, 0.0043737, 0.0045491, 0.0050336, 0.0057244, 0.0065365, 0.0074143, 0.0083245, 0.0092477, 0.0101727, 0.0110934, 0.0120064)
var <- seq(-2, 2.25, 0.25)

# plot data
plot.dat <- data.frame(b = b, se = se, var = var)
plot.dat$ub <- plot.dat$b + 1.96 * plot.dat$se
plot.dat$lb <- plot.dat$b - 1.96 * plot.dat$se

# make plot
plot3 <- ggplot() +
  geom_ribbon(data = plot.dat2, aes(x = var, ymin = lb, ymax = ub), alpha = 1, fill = "gray75") +
  geom_line(data = plot.dat2, aes(y = b, x = var), color = "black") +
  geom_hline(yintercept = 0, color = "black", lwd = 0.5) +
  plot_theme(size = 14) +
  ylab("Marginal Effect of Effectiveness") +
  xlab(NULL) +
  scale_x_continuous(limits = c(-2, 2.25), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-0.1, 0.12), expand = c(0, 0)) +
  theme(plot.margin = unit(c(1, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        plot.title = element_text(hjust = 0, margin = margin(b = 20)))

# make plot
plot4 <- ggplot(data = dat, aes(distance)) +
  geom_density(color = "black", fill = "gray75", alpha = 1) +
  plot_theme(size = 14) +
  xlab("Distance") +
  ylab("Density") +
  scale_x_continuous(limits = c(-2, 2.25), expand = c(0, 0)) +
  scale_y_continuous(limits = c(0, 0.6), breaks = c(0, 0.6), expand = c(0, 0)) +
  theme(plot.margin = unit(c(0.25, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 33, b = 0, l = 0)))

# model with expenditures, support, and influence
b <- c(-0.0271943, -0.0218375, -0.0163272, -0.0106632, -0.004845, 0.0011275, 0.0072543, 0.0135353, 0.0199702, 0.0265586, 0.0332998, 0.040193, 0.0472372, 0.0544314, 0.0617741, 0.0692638, 0.0768988, 0.0846772, 0.0925969, 0.1006555, 0.1088506, 0.1171793)
se <- c(0.0144925, 0.0135396, 0.0125886, 0.0116404, 0.0106972, 0.0097625, 0.0088421, 0.0079453, 0.0070872, 0.0062914, 0.0055951, 0.0050523, 0.0047316, 0.0046956, 0.0049677, 0.0055187, 0.0062908, 0.0072272, 0.0082847, 0.0094339, 0.010655, 0.011935)
var <- seq(-3, 2.25, 0.25)

# plot data
plot.dat <- data.frame(b = b, se = se, var = var)
plot.dat$ub <- plot.dat$b + 1.96 * plot.dat$se
plot.dat$lb <- plot.dat$b - 1.96 * plot.dat$se

# make plot
plot5 <- ggplot() +
  geom_ribbon(data = plot.dat, aes(x = var, ymin = lb, ymax = ub), alpha = 1, fill = "gray75") +
  geom_line(data = plot.dat, aes(y = b, x = var), color = "black") +
  geom_hline(yintercept = 0, color = "black", lwd = 0.5) +
  plot_theme(size = 14) +
  ylab("Marginal Effect of Expenditures") +
  xlab(NULL) +
  scale_x_continuous(limits = c(-3, 2.25), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-0.08, 0.16), expand = c(0, 0), breaks = c(-0.08, -0.04, 0, 0.04, 0.08, 0.12, 0.16)) +
  theme(plot.margin = unit(c(0.25, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        plot.title = element_text(hjust = 0, margin = margin(b = 20)))

# model with expenditures, distance, and influence
b <- c(0.078452, 0.073497, 0.068603, 0.0637704, 0.0589998, 0.0542914, 0.0496456, 0.0450628, 0.0405432, 0.036087, 0.0316946, 0.0273662, 0.0231018, 0.0189016, 0.0147658, 0.0106944, 0.0066874, 0.002745)
se <- c(0.0106744, 0.009396, 0.0081847, 0.0070605, 0.0060554, 0.0052195, 0.0046245, 0.0043483, 0.0044305, 0.0048347, 0.0054749, 0.0062665, 0.0071485, 0.0080824, 0.0090444, 0.0100198, 0.0109993, 0.0119771)
var <- seq(-2, 2.25, 0.25)

# plot data
plot.dat <- data.frame(b = b, se = se, var = var)
plot.dat$ub <- plot.dat$b + 1.96 * plot.dat$se
plot.dat$lb <- plot.dat$b - 1.96 * plot.dat$se

# make plot
plot6 <- ggplot() +
  geom_ribbon(data = plot.dat, aes(x = var, ymin = lb, ymax = ub), alpha = 1, fill = "gray75") +
  geom_line(data = plot.dat, aes(y = b, x = var), color = "black") +
  geom_hline(yintercept = 0, color = "black", lwd = 0.5) +
  plot_theme(size = 14) +
  ylab("Marginal Effect of Expenditures") +
  xlab(NULL) +
  scale_x_continuous(limits = c(-2, 2.25), expand = c(0, 0)) +
  scale_y_continuous(limits = c(-0.04, 0.12), expand = c(0, 0)) +
  theme(plot.margin = unit(c(0.25, 1, 1, 1), "cm"),
        axis.title.y = element_text(margin = margin(t = 0, r = 20, b = 0, l = 0)),
        plot.title = element_text(hjust = 0, margin = margin(b = 20)))

# save plot
pdf("Fjelstul_Carrubba_APSR_figure3.pdf", width = 9, height = 11)
grid.arrange(plot1, plot3, plot5, plot6, plot2, plot4, ncol = 2, nrow = 3, widths = c(4, 4), heights = c(4, 4, 2))
dev.off()

#################################################
# end R script
#################################################
